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Abstract. We overview the EDABI method developed recently and combining the 
exact diagonalization and ab initio aspects of electron states in correlated systems and 
apply it to nanoscopic systems. In particular, we discuss the localization-delocalization 
transition for the electrons that corresponds to the Mott-Hubbard transition in 
bulk systems. We show, that, the statistical distribution function for electrons in 
a nanochain evolves from the Fermi-Dirac-like to the Luttinger-liquid-like with the 
increasing interatomic distance. The concept of Hubbard subbands is introduced to 
nanoclusters, and corresponds to the HOMO-LUMO splitting in the molecular and 
organic solid states. Also, the nanochains exhibit magnetic splitting (Slater- like), even 
without the symmetry breaking, since the spin-spin correlations extend over the whole 
system. Thus, the correlated nanoscopic systems exhibit unique and universal features, 
which differ from those of molecular and infinite systems. These features define unique 
properties reflecting "the Mott physics" on the nanoscale. We also employ the EDABI 
method to the transport properties in nanoscopic systems. For example, we show 
that the particle-hole symmetry is broken when the tunneling conduction through H2 
molecule is calculated. 



PACS numbers: 61.64,+w, 36.40.-c, 73.63.-b 



Submitted to: J. Phys.: Condens. Matter 



The combined exact diagonalization- ab initio approach 



2 



1. Introduction 

The studies of nanosystems are becoming increasingly important in view of their 
application in quantum nanoelectronics and related fields of research. Of particular 
importance are their quantum electronic properties, since they determine their behavior 
as concrete devices: quantum nanowire connectors and semiconducting elements, 
single-electron transistors, spin valves etc. Under these circumstances, solid-state and 
molecular nanophysics is developing very rapidly to provide proper quantitative (and 
qualitative) characteristics of their static, transport, and optical properties [I]. 

Independently of the applications, the research in nanophysics is important for 
its own fundamental sake. Namely, the nanosystems are finite systems and therefore, 
most of the limiting situations considered in the condensed matter physics and involving 
the limit for the number of atoms N — > oo, is simply inapplicable. Furthermore, the 
role of boundary conditions is very nontrivial, as they should reflect the system actual 
configuration. Also, and probably most importantly, the question of electron-electron 
interactions, particulary for extended states, should be taken into account on equal 
footing with the single-electron aspects of their quantum states, since the screening 
processes are very often highly ineffective. Nevertheless, given the circumstances, the 
role of physics is to single out universal properties of the systems such as nanowires, 
clusters, quantum dots, etc. Such research involves also determining the conditions, 
under which the bulk-solid concepts are applicable, since then the description can be 
simplified remarkably. 

In this article we present our recent and earlier results concerning the studies 
of nanosystems containing as an intrinsic property the electron correlation effects 
induced by the Coulomb interactions, within the devised EDABI method combining 
both exact diagonalization and ab initio aspects of electronic states, as well as their 
transport properties [2]. This approach allows for discussing the evolution of the 
physical properties of e.g. nanowires or clusters, both as a function of their interatomic 
distance, as well as to determine the values of microscopic single-particle and interaction 
parameters. The visualization of the properies as a function of interatomic distance 
is particularly important for nanometer size systems, as they are studied customarily 
by placing them on a substrate with the lattice parameter, which differs from their 
equilibrium interatomic distance (sometimes the substrate even stabilizes them). 

Our most interesting results can be summarized as follows. First, we show how 
the system properties evolve from the Fermi-liquid-like to the atomic-like states for 
nanowires containing up to N = 16 simple atoms, passing through the Luttinger-liquid- 
like state, with the increasing interatomic distance R. This evolution is determined 
by calculating directly the system statistical distribution function n ka of electrons 
and their dynamical spectral function. Finite-size scaling properties are introduced to 
determine the Mott critical interatomic distance ( the Mott criterion ) for the transition 
from the delocalized to the localized states. Secondly, we show that the electronic 
states in nanosystems of iV ~ 10 atoms exhibit a magnetic splitting reminiscent of 
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the Slater splitting in antiferromagnetic metallic systems This type of symmetry 
change is associated with the breakdown of the discrete translational symmetry, as the 
antiferromagnetism sets in. In the case of nanosystems such symmetry breakdown is not 
required. It turns out, that the splitting appears if the spin-spin correlation length is of 
the system size. The splitting should be detected e.g. in nanowires containing strongly 
correlated electrons in a half-filled valence band configuration. Thirdly, the question 
arises whether with the increasing interatomic distance one should not observe the 
Hubbard split-band structure of nanowire and molecular (e.g. or Li n) clusters? We 
show that nanocluster levels group into multiple Hubbard subbands (in bulk systems 
such a grouping is termed as HOMO and LUMO structures). Finally, in the case 
of a quantum dot composed of e.g. H 2 molecule attached to the semimacroscopic 
electrodes we show, that the electron-hole symmetry in the tunneling transport through 
the molecule is not preserved and is due to the difference in electronic binding energy of 
H2 and i?2~ states. Such result is not obtained if a parameterized model of a quantum 
dot properties is used jlj. 

The structure of this paper is as follows. In the next Section we overview the EDABI 
method [SI IH El El Ej , as well as discuss some of the many-electron general properties 
adopted to the analysis carried out in the next Sections. We also compare the method 
with the configuration interaction (CI) approach used in quantum chemistry. In Section 
3 we discuss atomic systems and Hn nanoclusters, whereas in Section 4 the results 
for nanochains containing up to iV = 16 atoms are elaborated in detail. In Section 5 
we employ the EDABI method to calculate the tunneling conductivity through the H2 
molecule and the Drude weight for the nanochains. We also discuss the role of boundary 
conditions there. 

The present method and the research grew out of our earlier work on the 
thermodynamics of the Mott-Hubbard transition in correlated systems jH]. There, the 
principal question was if the metallic and magnetic insulating states can be regarded 
as separates phases in the thermodynamic sense? The affirmative answer to the above 
question provided a partial answer to the Sir Nevill Mott question: What is a metal? In 
this article the corresponding principal question: how small a piece of metal can be? In 
other words, can a nanochain composed of e.g. = 10 Cu atoms be regarded already 
as a metallic system and in what sense? The answer is 'y es \ but above a critical value 
of voltage applied to the chain and for not too large interatomic distance. 

2. Exact diagonalization combined with an ab initio approach 

2.1. General features 

The second-quantization language is used when we have an interaction between the 
quantum physical fields representing the classical particles composing the system. The 
approach is usually formulated in the occupation-number representation, expressing the 
possible occupations of a given {complete in quantum mechanical sense) set of single- 



The combined exact diagonalization- ab initio approach 



4 



particle states, between which there are transitions induced by their mutual interaction. 
Explicitly, the single-particle basis {<3>j(r)} defines the field operator ^(r), in terms of 
which the many-particle Hamiltonian (or Lagrangian) is defined. In the nonrelativistic 
case, the Hamiltonian is defined as [5] 

H = J rf 3 #(r)ffi(r)*(r) + i J d 3 rd 3 r'& {r)&(r')V{r - r')£(r')#(r) 
= f + V (1) 

with 

^(r) = ^$ J (r)a J . (2) 

i 

H\(r) represents Hamiltonian for a single particle in the assembly of N indistinguishable 
particles, V(r — r') is the interaction between a single pair of particles, aj is the anni- 
hilation operator of the particle in the single particle state $i(r). One should underline 
that the basis {$j(r)} should be complete, but otherwise arbitrary (it does not have to 
be orthogonal [10). 

The complete set of {<3>j(r)} in (j2J) is in practical calculations not infinite, so we 
have to resort to a finite subset of [$i(r)] when performing explicit computations, i.e. 
constructing many-particle models. In this manner, we make the basis incomplete in 
the quantum-mechanical sense, even though the model may be well justified on physical 
grounds as it may contain principal dynamic processes involved. In the first step, we 
diagonalize H in the Fock space for given trial basis [$j(r)]. Then, as a next step, 
one can optimize the finite (incomplete) subset taken by e.g. minimizing the ground 
state energy E G = ($ \H\Q> ) (where |<3> ) is the ground state wave function in the 
Fock-space representation), with respect to the selected subset [$j(r)]. Such procedure 
is highly nontrivial, since the determination of Eq = ^{^(r), V$j(r)} requires first 
the diagonalization of the parameterized Hamiltonian (|TJ in the Fock space pi] (with 
the microscopic parameters containing both [$i(r)] and [V$i(r)]) and, only after that, 
setting up an effective wave equation for each $j(r), with Eq treated as a functional of 
the trial basis [$*(r)]. The resultant renormalized or self-adjusted wave equation (SWE) 
should have a universal meaning to the same degree, as has the starting Hamiltonian (0) 
in that incomplete basis. We call this method of approach EDABI ( a combined Exact 
Diagonalization - AB Initio approach ). In this manner, only then the approach to 
the interacting system can be regarded as completed, particularly in the situation when 
the interparticle interaction cannot be regarded as weak, e.g. for correlated systems. 
Obviously, when the correlations are weak, the approach reduces to the Hartree-Fock 
approximation. This approach has been implemented so far to some exactly soluble 
models of correlated electrons and to nanochains |T2J EH E] containing up to N = 16 
atoms using as a trial basis, with adjustable Slater or Gaussian orbitals j3j. Here we 
present a rather extensive analysis of our method, applicable to both fermion and boson 
systems, as well as construct explicitly the multiparticle wave function is the simplest 
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atomic situations. The last aspect of the work may be applied to both atomic and molec- 
ular systems providing e.g. a systematic approximation scheme in the quantum-chemical 
calculations. In particular, the explicit many-particle wave-function construction allows 
for a comparison with the method of multiconfigurational interaction (MCI) approach 
[TK] utilized in quantum chemistry. Simply put, we develop a relatively straightforward 
(but not simple!) workable scheme, which is applicable to both many- electron atoms, 
molecules and molecular ions, as well as to clusters and nanoscopic systems. 

One can summarize important features of our approach as follows. First, the inter- 
action and the single particle terms in the many-particle Hamiltonian are treated on an 
equal footing. Second, the single-particle wave function appears in a nonexplicit form 
in the expression for the microscopic parameters and is determined explicitly from the 
variational principle for the ground-state energy Eg = (T) + (V) = Ec{&i(r)}, which 
leads to the self-adjusted wave equation [IB]. Third, the explicit construction of the 
multiparticle wave function ^(rj, ...,r N ) in terms of the field operators [T7j provides a 
way of systematizing the approximation schemes. Fourth, one does not count twice the 
interaction, which should be taken with care in e.g. in LDA+U method. All of these 
features are particularly important for the systems, for which the interaction cannot 
be treated as a perturbation. Among such systems are cluster and various correlated 
Fermi and Bose systems, particularly of low dimensionality. Probably the most impor- 
tant formal feature of the present approach is that that by reversing the order of solving 
the problem (treating first the interaction in the Fock space and only then determining 
the single particle wave function in the Hilbert space), we complete the treatment of 
strongly correlated systems. Unfortunately, our approach is executable so far only in 
limited number of situations. As a side result we obtain also the values of microscopic 
parameters for the parametrized models [18] . In effect, the physical properties can be 
analyzed as a function of interatomic distance, not only as a function of model param- 
eters, and thus provide us with the global minimum for system at hand and for given 
interatomic distance. 

The self-adjusted wave equation (SWE), as we shall see, is of nonlocal and nonlinear 
nature. Hence, it is very difficult to solve it directly. Nonetheless, the main purpose of 
the present paper is to present the solution in the closed variational form and illustrate 
its character in simple situations, ranging from the atomic physics to nanophysics. 

2.2. Renormalized (self-adjusted) wave equation (SWE) 

We start with Hamiltonian (JIJ and write down in the explicit-spin basis, in which the 
spin is regarded as an additional coordinate, i.e. in the form 




+ / /^'i^i 1 ( r i)^! a ( r a)V'(ri-r a )* asi (r a )i <71 (r 1 ). (3) 
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We define the spin-dependent field operator as 



M 



^<x( r ) = ^2w t (r)xaaia, (4) 
1=1 

where {wi(r)} is a complete single-particle basis with the set of quantum numbers 
denoted by i. Note that we regard the Hamiltonians Hi(r) an d V(ri — r 2 ) as spin 
independent (it is straightforward to generalize the formalism to the case with spin- 
dependent Hamiltonian, i.e. when magnetic field or spin-orbit interaction are included). 
Substituting (j3J) into Q we obtain the usual form of the Hamiltonian 



H - ^2tija\ a a jcT + - Vijkl<4a 1 a ja 2 a lv2 a kari, (5) 

ijcr ijkl(j\(J2 

with the microscopic parameters defined by 

Uj = (wi\Hi\Wj) = J dhwti^H^Wjir), (6) 

and 

Vijki = (wiWj\V\w k wi) (7) 

ee J dh 1 d 2 d^r 2 wl{r 1 )w](Y 2 )V{r 1 -r 2 )w k (Y 1 )wi(r 2 ) 

In the standard form © of the many-particle Hamiltonian the single- and many- particle 
aspects of the problem are separated in the sense that the calculation of the hopping 
parameters^.,- and their interaction correspondants V^ki, both containing the single- 
particle wave- functions {wi(r)}, is separated from the diagonalization procedure of the 
Hamiltonian in the Fock space ^Hj- The latter procedure is dependent only on the 
nature of the commutation relation between the annihilation (a ia ) and creation 
operators. Thus this two-step procedure can be seen explicitly when we calculate the 
system ground-state energy 

E G EE (H) =^2^(4^) + - ^ V m{ a \aA<r2 a ^ a k*l) >> ( 8 ) 
iju ijkla\c>2 

where the averaging (...) takes place over all accessible occupancies of given single 
particle states \iai), | jcx 2 ) , \ko~i), and |Zcr 2 ). Obviously if we want to consider all 
occupancies (a grand canonical ensemble), we diagonalize H — fiN, where N is the total 
number of particles and only a posteriori impose the conditions that iV = YliA n io)i 
with n ia ee a\ a a ia 

So far, the approach is standard [HI ED] • We have proposed [21 H2 EH EI to close 
the solution (i.e. the complete calculation of e.g. Eg) with the determination of the 
single-particle basis {wi(r)} by treating the expression (JBJ as a functional of the set of 
functions {u>;(r)} and their gradients. In such situation the renormalized (self- adjusted) 
wave equation is determined from the Euler equation for the functional 

E{ Wi (r)} ee E G {w t (v)} -^N-J2 *ii ( / d 3 rwt(r) Wj (r) - <Sy J , (9) 

i>j \ J ' 
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where 



A" is the number of particles in the system, and \j are the Lagrange multipliers, when 
the single-particle basis is nonorthonormal. 

The general form of this equation in the stationary case and in the grand canonical- 
ensemble formalism is 

We will make a fundamental postulate concerning this equation: ^4s this equation 
does not contain explicitly the (anti) commutation relations between the creation and 
annihilation operators, it is equally valid for both fermions and bosons and determines a 
rigorous, within the class of states included in the definition o/ 1 l r (r) ; wave equation for 
a single particle wave function in the ground state, in the millieu of remaining (N — 1) 
particles. Additionally, as is implicit in the treatment above, we have defined one global 
spin-quantization axis for all single particle states used to define ^ CT ( r )- I n some spin 
noncolinear systems this is insufficient and will require a more refined treatment. Also, 
if we use the particle conserving approach to calculate E G , then we put /x = in (fTUj) . 
Likewise, for the orthonormal basis {uii(r)} used to define ^(r) in (jlj), one puts Ajj = 0. 
In the latter case the system (fTtH) represents a set of Euler equations for renormalized 
Wannier functions. In what follows we discuss examples of application of this equation 
to fermion systems (it can be applied to Bose systems in the same manner). But first, 
we define the renormalized many particle wave function as complementing the above 
renormalized single particle states {uii(r)} and then discuss the differences with the MCI 
approach. 

2.3. Multiparticle wave function from the second-quantization approach 

The general N-particle state |$jv) in the Fock space can be defined through the 
corresponding A^-particle wave function \& a (ri, . . . , tn) in the Hilbert space in the 
following manner [01 El 

\$ N ) = ^=J d 3 ri...r N * iV (r 1 ,...,r N )it( ri )... 1 |rt( rN )|o) ) (n) 

where |0) is the vacuum state. One can reverse this relation and a simple algebra yields 
the following expression for the wave function ^(ri, . . . , rjy) in terms of |$jv) 

^(r 1; . . . , r N ) = _L<0|*(rO . . . *(r N )|$*>. (12) 

In other words, to obtain the wave function in the coordinate representation, we not 
only annihilate N-particles from the state |$jv), but also project out thus obtained 
result onto the Fock vacuum state and normalize it by the factor (AH) -1 / 2 . Usually, 
the formula (fT2"|) is not used as we proceed from first to second quantization. Now, 
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Figure 1. Flow-chart describing the scheme of the EDABI method. For details see 
main text. When selecting the Gaussian starting single-particle set, the topmost block 
should be disregarded. 



the crucial point is based on the observation that if we substitute in the field operator 
\l/(r) the renormalized wave functions obtained from Eq. (fTU|) . then we should obtain the 
renormalized field operator and as a consequence, the renormalized multiparticle wave 
function \I/ a (ri, . . . , tn) from relation ()12j) . This last step of inserting renormalized 
field operator completes the procedure of a formal treatment of many-particle system, 
which avoids writing down explicitly the N-particle Schrodinger equation. The whole 
approach is schematically represented in Fig. ^ This scheme provides an exact 
renormalized single-particle wave function from Eq. and the exact iV-particle wave 
function provided we can diagonalize the second-quantized model Hamiltonian (jSJ) for 
the problem at hand. We shall see next, that we can approach the true solution also 
step by step. 

2.4- Finite-basis approximation for the field operator: Difference with the 
multiconfiguration interaction (MCI) approach 

The field operator \l/(r) defined in terms of the sum over a complete basis {wi(r)} 
contains an infinite number of single-particle states. We assume that, in general, we 
represent the field operator by M wave functions {u>j(r)}. Explicitly, 

co M 

^( r ) =^2 w i( r ) a i-^2 w i( r ) a i: ( 13 ) 

i=l i=l 

with i representing a complete set of quantum numbers and M being a finite number. 
This approximation represents one of the most fundamental features of constructing 
theoretical models. The neglected states usually represent highly excited (and thus 
negligible) states of the system. We can then write the approximate N-particle wave 
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function (N < M) in the following manner 

1 M 

* Q (ri, . . .r n ) = — = ^2 (°\ a iN ■ ■ •Oiil^Jv)^ il (ri) . . .w ijv (r N ). (14) 

V ' ii,...,ijv=l 

Recognizing that within the occupation-number space spanned on the states 
{\ik)}k=i...M, we have the N-particle state in the Fock space of the form 

1 - 

l^) = ^ft E C«4v4l )' ( 15 ) 



v'N\ . 

JlvJiV = l 

where Cj lt „j N represents the coefficients of the expansion to be determined from a 
diagonalization procedure. Substituting (JT5j) to (jHJ) we obtain 

^ a (n,...,r N ) = (16) 

M M 



1 



V (0|a n ...a iJv aj i ...a J t iv |0)C jl ... ijv w il (ri)...w ijv (r N ). 



AH ^ 

h,...,ijv=1 ji,...,ijv=i 

The expression provides AH nonzero terms each equal to (— l) p , where P represents 
the sign of the permutation of quantum numbers (ji-.-jjv) with respect to selected 
collection . . . z'jv)- In other words, we can write that 

1 M 

tf a (ri,...,r N ) = — Ci 1 ...iw(A^)[«'< 1 (ri)...'u; iAr (r N )]. (17) 

ii,...,ijv=l 

We have the same expansion coefficients for both wave function in the fock space 
|<&jv) and that in Hilbert space\E' Q ,(r 1 , . . . , r N )! Therefore, the above expression 
represents the multiconfigurational-interaction wave function of iV particles distributed 
among M states with the corresponding weights C^...^ for each configuration, and 
(A, S) represents respectively the antisymmetrization (Slater determinant) or the 
symmetrization (simple product u^ri) . . . Wi N (r-^)) for the fermions and bosons, 
respectively. Whereas MCI used in quantum chemistry [20] bases on variational 
optimizations of both the coefficients C^...,^ and the basis {u>;(r)}, here the coefficients 
C are determined from diagonalization in the Fock space, spanned on M states in the 
Hilbert space and determined from (fTt)|) . The presence of SWE (|T0|) thus supplements 
the usual MCI approach. 

Summarizing, the differences between the EDABI and the MCI methods, both of 
which belong to the class of multi-determinant expansion of N-particle wave function, 
are threefold: 

i) Historical. MCI evolved from variational methods of quantum physics and 
chemistry to include the electronic correlations and hence, to obtain lower value 
of Eq by starting from many-particle Schrodinger equation. EDABI represents a 
procedure of calculating single-particle wave function starting from parameterized 
models of strongly correlated electrons. 
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ii) Technical. In MCI, we optimize simultaneously the coefficient expressing the 
weights of different determinants (representing different micro-configurations), as 
well as the parameters of the trial single-particle basis. In EDABI, we diagonalize 
the Hamiltonian expressed in the Fock space (with the help of either analytic or 
numerical methods), combined with a simultaneous optimization of the orbital size 
in the resultant ground state. 

Hi) Essential. In the case of analytically soluble models, EDABI leads formally to 
the explicit form of the renormalized wave equation, which represents a nonlinear 
Schrddinger equation of nonlocal type. This circumstance opens up a new direction 
of studies in mathematical quantum physics. Additionally, it allows for a direct 
determination of dynamical correlation functions, transport properties, etc. in the 
convenient, second-quantization, language. 

2.5. An interpretation of the approach: rationale behind the self-adjusted wave function 
in the many -body system 

Usually, the choice of starting single-particle basis {wi(r)} is dictated by the physics 
of the system at hand. Since the creation and annihilation processes are characterized 
only by the quantum numbers i of those starting single-particle states, one can say 
that they represent a particle language characterizing transitions between those states. 
The ground state energy obtained from the diagonalization in the Fock space defines 
resultant single-particle states, which can be called the self-adjusted states, after the 
optimization of the single-particle wave function has been carried out via solving the 
self-adjusted wave equation (SWE)j2HI22| or its variational version. In other words, we 
allow the initial particle wave function Wi(r) to adjust to the correlated state. In such 
scheme the particle and the wave aspects of the single-particle states are intertwinned 
formally, illustrating among others the particle-wave complementarity, this time in a 
formal manner. For example, the Born probability density of finding a particle is taken 
here as ^] o .(^j T (r) x I' -(r))/A^. In what follows we essentially illustrate the method on 
concrete examples and apply it to nano systems. 

One may ask the basic question: Why to revert the usual sequence of solving the 
single-particle wave equation first, and only then constructing the field operators in the 
second-quantization Hamiltonian (or Lagrangian), by solving the many-particle Hamil- 
tonian first and only then readjusting the orbitals? The reason for this is as follows. As 
said above, in many cases (e.g. for correlated and/or low-dimensional fermionic systems) 
we encounter the situation when the interaction cannot be regarded as a perturbation 
and therefore should be treated at least on equal footing with the single-particle aspect 
of the problem. This is because, in general, the interaction may change the class of the 
stationary-state wave function. Such a situation is beautifully illustrated on example of 
metal-insulator phase transition of the Mott-Hubbard type, at which the metallic state 
represented by the Bloch-type wave functions switches to the localized (Wannier-type) 
states even though those representations are regarded as equivalent from a single-particle 
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point of view. In other words, the interaction determines the particle wave function. 
Furthermore, the electron as a separate entity (lepton) is preserved even in the highly 
correlated milieu of other particles and therefore, constructing the self-consistent wave 
equation (SWE) has a sense as it provides its wave function adjusted to the environment. 



2.6. Single-particle basis selection and the particle- density space profiles 

As we have already mentioned, the selected single-particle basis Wi(r) = {$j(r)} is 
determined from the variational principle for ^{^(r), V$i(r)}, and it should satisfy 
the self-adjusted wave equation. Here we solve this equation variationally and take a 
trial basis {$j(a;r)} dependent on a finite number of parameters a = {a p } p =i^K- 
Moreover, if the basis is orthonormal, then the equation can be simplified then to the 
form 

5(E G - uN) , s 

V G , ; = 0. 18 
oa 

In most applications we select adjustable atomic orbitals as the basic functions 
r)}, which need to be orthonormalized with a special kind of orthonormalization 
procedure. Additionally, if the number of particles N is conserved by the Hamiltonian 
H then the term fiN in (118)1 is absent, and the equation reduces to the ordinary energy 
minimization E G = (H) with respect to the variation parameters {a p }. 

2. 7. A remainder: Wannier basis for an extended periodic system 

Let us consider the multi-particle systems for a single-band case for atoms located at 
positions {R«} and the basic wave functions r)} of the form r — Hi)}. The 

orthonormalization procedure for periodic structures can be obtained by starting with 
the expression for the Bloch function 

N 

$ q («;r) = iV q ^ e ^$ J («;r), (19) 
3=1 

where iV q is a normalization coefficient. The normalization condition then takes form 
($ q («;r)|$ q ,(a;r)) = N^N^e^^^-^ia^^T)) 

where is an overlap integral. The last sum does not depend on the relative distance, 
hence 

($ q («;r)|$ q ,(«;r)> = N^NS^ ]T e^ R - R ^ = <W, (21) 

1 

and thus 

I A/q| = (Nj2 ei(l<ni ~ nj)s i3)~ 1/2 - ( 22 ) 
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As we can see, the above wave functions {<E> q (a; r)} are the Bloch functions, from 
which we can construct the Wannier functions {wi(a;r)} in the usual manner, i.e. 



^(«;r) = ^E q e- iq ' Ri $ q (a;r) 

(23) 



^E q e- iq - Ri Ar q E,e^^(a;r) 



where 



= E,^(«;r)(^E q ^e-^( R -^)) 

i p -iq-(Ki-nj) 

Pu = — Y , (24) 
3 V / E m e- iq ' (R "- Rm) ^ 

The orthonormal Wannier set {wi(a;r)} obtained with the help of the above 
procedure is thus a set of linear combinations of the originally non-orthonormal atomic 
orbitals. Moreover, the coefficients {Pij} depend only on R; — Ft, and satisfy the relation 

2. 8. Wannier functions for finite systems 

The above procedure can be applied with modifications to finite cluster systems. We 
start from the decomposition 

M 

wi(a;r) = ^2l3ij^ j (a;T), (25) 
j'=i 

with the normalization condition in the form 

(u> z (a;r)K(a;r)> = £\ fc (3^,^ (a; r)|$ fc (a; r)> 

= YjjkPljPl'kSkj = Sll>- 

In the matrix language, this condition can be rewritten as 

(3S(3 + = 1, (27) 
or, equivalently as 

= s- 1 . (28) 

We choose the (3 matrix in the form (3 = (3 + ; in effect this choice leads to the relation 

(3 = S~ 1/2 . (29) 

The above method is known as the Lowdin method of determining the orthonormal 
basis, and for the overlap integral matrix becomes unity if elements of the system are 
separated from each other at large distances. So, the (3 matrix then can be written down 

as 

00 / i\n An n co n 

/'= 1 + E 1 ^n( 2 *- 1 ) = 1 + E A "n(^- 1 ). < 3 °) 

n=l k=l 7i=l k=l 
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where A = 5 — 1. However, the above series does not converge in the tight binding 
approach and it has to be modified to the form 

and thus 

oo n 1 

(3 = (1 + CT 1/2 [1 + J2 A c - !)]» ( 32 ) 

n=l k=l 

where A c = S/(l + C) — 1. The parameter C allows us to manipulate convergence of 
the series, for it appears in the convergence condition 

l|Ac|| < 1, (33) 

that needs to be satisfied for any kind of the || • || metrics. Therefore, we choose the 
metrics 

|| Aclloo = max{|^- - 5 t A} = max{-^— ; max ^Jf^ }, (34) 
11 ciioo „ \\ 1 + c *j\s X l + C 1 + C 1 V ; 

As we can see, the metrics satisfies the condition when the parameter C is equal to 

C = max \ Sij\. (35) 

The method, we have just shown, can be applied to both single-band and to multiple- 
band systems. The method can also be applied to disordered systems. In what follows 
we will select as a starting atomic basis {$j(r)} either Slater or Gaussian (STO-3G) 
basis when constructing the Wannier function, and subsequently, optimize their size in 
the correlated state. 

2. 9. Particle density profiles in space 

In the previous Section we dealt with the N-particle wave function ^(ji, ...,r N ). Here 
we show how it can be applied to the evaluation of particle density n(r). We start with 
the usual definition of the probability density for a single particle: 

p(r N ) = y'd 3 n...r Ar _i|*(ri,...,rj V )| 2 . (36) 

We utilize now the expression (4) for the wave function we obtained in Section II. In 
effect, 

|^(r 1 ,...,r w )| 2 = ^($ J v|$ + (r JV )...$ + (r 1 )|0)(0|§(r 1 )...§(r w )|<I> J v). (37) 



We have the particle-number conservation and hence we can rewrite (37) as 

1 

W. 



^f(r 1 ,...,r N )\ 2 = ■4<$Jv|* + (r JV )...$ + (ri)$(ri)...$(r JV )|$ JV ). (38) 
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This is because we can insert J2 N \N)(N\ instead of |0)(0|. The expression for the field 
operator \I/(r) , defined with the help of the orthonormal basis {u>j(r)} in (4), leads to 
the following relation 

$ + (r)$(r)= <( r )x^ 2 WXa 2 < CT1 <W 2 , (39) 

il<7li2<T2 

and thus 

¥ + (r)$(r) = £ wUv)w i2 {v)at lC a i2a . (40) 

Therefore, we can apply the above relation to the expression for pir^). Namely, by 
noting that 

p(r N ) = J d 3 r 2 ...r N ^( J d 3 n|tt(n, r N )\ 2 ), (41) 
we can obtain the following series of helpful identities: 

f^l^n,...,^)! 2 

= i / rf 3 r 1 ($^|§+(r iV )...$+(r 1 )§(r 1 )...§(r i v)|$ J v) 

^Ei^O^'iKfaK^ (42) 
= ] i T ($ J v|$+(r JV )...§+(r 2 )E lCT n l(7 §(r 2 )...§K)|$ JV ) 
= ^($ JV |$+(r iV )...$+(r 2 )$(r 2 )...$(r iV )|$ Ar ). 

The expectation value of the particle-number operator ^2 ia n ia is equal to unity for the 
multi-particle state ^(r 2 )...^(rjv)|$jv) ~ |$jv=i)- Obviously, the subsequent procedure 
applied N — 2 times leads to the net result 

p(r N ) = ^($ N \* + (r N )V(r N )\$ N ). (43) 
This result can be understood if we introduce explicitly the particle-density operator 

n(r) = $ + (r)§(r), (44) 
and thus the particle density is 

n(r) = Np(r) = ($ N \* + (r)$(r)\$ N ). (45) 
Hence, we obtain the explicit expression for the density of particles using Eq.(4), 

i.e. 

n(r) = <( r ) w *2(r)(®N\at ia a V2a \$ N ), (46) 
or, equivalently 

n(r) = ^ kWI 2 ^) + ^ <(r> i2 (r) ^« CT a, 2(7 >, (47) 

where averages (...) are taken with the N-particle ground state \<&n)- 

The first of two terms represents the contribution of the particle-number operator to 
the total density of particles and appears in the Hartree-Fock approach as the only term. 
The second term can provide a significant contribution when the averages {af ia a i2a ) is 
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of the same magnitude as the first contribution. Explicitly, if the N-particle function 
has the form of a simple determinant 



1 



H...IJV 



HI; 



W; 



(48) 



^(TnJ ... w iN {r N ) 

then the occupancies (n») = 1 in (|l*7jl an (a\ \ a a>i 2 o) = for z'j 7^ z 2 . This is not the case 
when the multi-configurational form (|T2"j) of ^(i - !, . . . , tn) is taken, as we shall see in 
the following. In the remaining part of this paper we discuss the results obtained within 
this method of approach for various nanosystems. 



3. Atomic systems and nanoclusters 

In this Section we start with the simplest examples of lightest atoms and ions to illustrate 
the specific features of the EDABI method, as well as to provide an elementary example 
of the renormalized wave equation. In particular, we systematically enrich up the trial 
basis to build up the Fock space. 



3.1. A didactic example: He atom 

We start by selecting as {wi(r)} just two Is- like states for the He atom $ CT (r) = 
(a 3 /^) 1 / 2 exp(— ar)xai where a is the effective inverse radius of the states. In other 
words, the simplest trial field operator is of the form 

$(r) =$ T (r)a t + $1(^01, (49) 

where a a is the annihilation operator of particle in the state $o-(r). The Hamiltonian 
in the second quantization for this two-element basis has then the form 

H = e a (n T + nj_) + Un^n^ (50) 

where nj = atof, whereas 

e a ={$ a \H 1 \$ a ), (51) 

and 

U = (^ w \V\^ a ) (52) 

are the matrix elements of the single-particle part defined as 

^ 2 ^2 ^ 2 ^2 2e2 2e2 «•«■ ^2 ^2 4 4 „m 
Hi = ~^V 2 ~ 77- V2 = -V 2 - V 2 . 53 

2m 2m kqTi kqT2 r\ r 2 

and of the Coulomb interaction 

V = -r^ r 1— ^— r. ( 54 ) 

Ko|ri-r 2 | |ri-r 2 | 
with the corresponding definitions in atomic units after the second equality sign. The 
only eigenvalue of (J5UJ) is obtained for the state ajajjO > and is E = 2t a + U. This 
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total energy is then minimized with respect to a to obtain the well-known variational 
estimate of both a and the ground state energy E G , as discussed before 0. However, we 
may look at the problem differently. As the approximate field operator can be defined 
in an arbitrary basis, we may regard the eigenvalue E as a functional of $o-(r), since 
the functions are under the integral expressions. Therefore, the true wave function is 
obtained from the Euler equation for the functional under the proviso that the wave 
function is normalized. This means that we minimize the functional 

£{$ CT (r)} = J2 f ^Cf^iWiW (55) 
+\J2[ dWr'|<Mr)|V 12 (r-rO|<Mr)| 2 . 

a 

In effect, the Euler equation take the form of he unrestricted Hartree-Fock equations 
for $ CT (r) 

(V - —) $ CT (r) + $ a (r) / d 2 r' f J MOI 2 = A$,(r). (56) 
V K o r r J J K \r-r'\ 

Thus we can see that taking in the simplest case just two spin orbitals we obtain either 
well-known variational estimate [21] for a and Eq for He atom: a = 27/(16ao) and 
Eq = — 5.695-Ry, where Oq ~ 0.53A is the Is Bohr orbit radius. 

Obviously, the proposed expression (J49j) of the field operator is the simplest one, 
though it leads to nontrivial results even though the trial atomic basis {$ CT (r)} is 
far from being complete in the quantum-mechanical sense. However, we can improve 
systematically on the basis by selecting a reacher basis than that in ([49)1 . The further 
step in this direction is discussed next. 



3.2. The basis enrichment for the lightest atoms and ions: He, Li 

We can expand the field operator in the basis involving the higher order irreducible 
representations of the rotation group with n = 2, which in the variational scheme 
involve including, apart from the ^i s (r) orbital, also the higher *f?2s( T ) an d \&2pm( r )> 
with m = ±1,0 (i.e. the next shell); all of them involving the adjustment of the 
corresponding orbital characteristics i = Is, 2s and 2pm. The field operator is then 

+i 



m=—l 



^2wi(r)xiaaia, (57) 



where u>j(r) are orthogonalized orbitals obtained from the nonorthogonal atomic 
basis {\I/j(r)} in a standard manner. The Fock space spanned on 2 + 2 + 6 = 10 trial 
spin orbitals contains D = dimensions, where M = 5 now and N e = 2,3 is the 
number of electrons for He and Li, respectively. This means that D = 45 and 120 in those 
two cases and we have to diagonalize the Hamiltonian matrices of that size to determine 
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Table 1. Optimized Bohr-orbit radii at — a~ of Is, 2s, and 2p orbits (in units of do), 
the overlap S between renormalized Is and 2s states, and the ground state energy for 
the lightest atoms and ions (five Slater orbitals taken). 





an 


a 2s 


&2p 


S 


E G (Ry) 


H 


1 


2 


2 





-1 


H~ 


0.9696 


1.6485 


1.017 


-0.1 


-1.0487 


He 


0.4274 


0.5731 


0.4068 


-0.272 


-5.79404 


He- 


1.831 


1.1416 


0.4354 


-0.781 


-5.10058 


Li 


0.3725 


1.066 


0.2521 


0.15 


-14.8334 


Be+ 


0.2708 


0.683 


0.1829 


0.109 


-28.5286 



the ground and the lowest excited states. One should note that we construct and 
subsequently diagonalize the < i\H\j > matrix in the Fock space for (fixed) parameters 
e a ,tij, and Vki- After the diagonalization has been carried out, we readjust the wave 
function and start the whole procedure again until the absolute minimum is reached (cf. 
Fig. HJ). 

By diagonalizing the corresponding Hamiltonian matrices and subsequently, 
minimizing the lowest eigenvalue with respect to the parameters a« - the inverse radial 
extensions of the corresponding wave functions, we obtain the results presented in Table 
^ (the values a2 pm are all equal within the numerical accuracy ~ 10~ 6 ). For example, 
the ground state energy of He is Eq = —5.794 Ry, which is close to the accepted "exact" 
value [21] —5.8074, given the simplicity of our approach. Further improvement is feasible 
by either including the n = 3 states or by resorting to the Gaussian basis; these are 
not performed in this Part (see next Sections). Instead, we elaborate on two additional 
features. 

First, we can represent the ground-state two-particle spin-singlet wave function for 
He atom taking \&(r|) in the form (|57j). which has the following form jH] 

K e >= (-0.799211a+ si a+ T + 0.411751a+ 4 a+ T - 0.411751a+ T a+ 4 (58) 
-0.135451a+ 4 a+ T + 0.0357708a+ 0j a+ 0? + 0.0357641a+ i; a+ _ 1T 
-0.0357641a+ 1T a+_ u )|0>, 

Similarly, the S z = +1/2 state for Li atom is of the form 

>= (0.997499a+ 4 a+ T a+ T - 0.0570249a+ T a+ 4 a+T ( 59 ) 
+0.0039591a+ t a+ oi a+ OT + 0.00395902o+ T a+ u a+ _ 1? 
-0.00395894a+ T a+ 1T a+_ u - 0.0237830+^0+^ 
-0.0237806a+ T o+ u o+_ 1T + 0.0237806a+ T a+ 1T o+_ u )|0 > . 

We see that the probability of encountering the configuration Is 2 in He is less than 
2/3, whereas the corresponding configuration ls 2 2s for Li almost coincides with that 
for the hydrogenic-like picture. The reason for the difference is that the overlap integral 
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Table 2. Microscopic parameters (in Ry) of the selected atoms and ions all quantities 
are calculated for the orthogonalized atomic states, t is the ls-2s hopping magnitude, 
Ui is the intraorbital Coulomb interaction (i=ls(l), 2s(2), m=0(3), and m=±l(p)), 
whereas and Jij are the interorbital Coulomb and exchange interaction parameters. 

t Ui U 2 U 3 Up K 12 K 13 K 23 J 12 J 13 J 23 

H~ 0.057 1.333 0.369 0.77 0.728 0.519 0.878 0.457 0.061 0.138 0.035 

He 1.186 3.278 1.086 1.924 1.821 1.527 2.192 1.289 0.212 0.348 0.115 

He~ -1.1414 1.232 0.764 1.798 1.701 0.929 1.421 1.041 0.269 0.28 0.102 

Li -0.654 3.267 0.533 3.105 2.938 0.749 3.021 0.743 0.06 0.606 0.014 

Be + -0.929 4.509 0.869 4.279 4.049 1.191 4.168 1.175 0.105 0.837 0.025 



between Is and 2s states S =< ls\2s > in the former case is large and the virtual 
transitions Is f± 2s do not involve a substantial change in of the Coulomb energy. Those 
wave functions can be used to evaluate any ground-state characteristic by calculating 
< ^gIOI^g > f° r O represented in the 2nd quantized form. For example, the atom 
dipole moment operator is d = e J <i 3 r^'"'"(r)x^(r), etc. 

The second feature is connected with determination of the microscopic parameters 
Vijki in our Hamiltonian, since their knowledge is crucial for atomic cluster calculations, 
as well as the determination of physical properties of extended systems as a function of 
the lattice parameter. Namely, we can rewrite the Hamiltonian (jSJ) for the case of single 
atom within the basis (|57|) in the form 

5 j 

H = eiUia- + ( a 2 CT a i<x + aL«2 CT ) + ^2 Uin^riii + - ^ K^mnj 

ia a i=l i^j 

~\Y1 J v ( Si ' S J ~ \ Uin n + J ii a lAl a ii a rt + %^at CT a iCT .(60) 

t is the hopping integral between Is and 2s states, Ui are the intraorbital Coulomb 
interactions, Kij are their interorbital correspondants, Vy is the so-called correlated 
hopping integral, and is the direct exchange integral, for states % and j = 1, . . . , 5. 
The principal parameters for the atoms and selected ions are provided in Table El 
We can draw the following interpretation from this analysis. The calculated energy 
difference AE for He between the ground state singlet and the first excited triplet 
is -2.3707 - (-5.794) ~ 3.423% ( the singlet Is t 2s I is still 1 Ry higher). The 
corresponding energy of the Coulomb interaction in Is 2 configuration is U\ = 3.278, the 
value comparable to AE. Additionally, the Coulomb interaction in Is f 2s J, state is 
~ l.bRy, a substantially lower value. The relative energetics tells us why we have a 
substantial admixture of the excited Is j 2s j state to the singlet Is 2 . In other words, 
a substantial Coulomb interaction ruins hydrogenic-like scheme, although the actual 
values could be improved further by enriching still the trial basis. 

One may ask how the renormalized wave equation would look in the present 
situation? The answer to this question is not brief already for the basis containing 
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2 4 6 8 10 

interatomic distance R/^ 



Figure 2. The level scheme of the H2 ground state and the lowest states as a 
function of the interatomic distance R. The hopping electron illustrates the relevance 
of H2 ionic configuration when measuring tunneling conductivity of H2 system. 

M = 5 starting states {iUj(r)}; we return to this question in a slightly simpler case of 
molecular states we consider next. 

3. 3. H 2 molecule and the H ^ and ions 

In this Subsection we consider H 2 and Li 2 molecules, and H 2 ion molecule. For the 
illustration of the method we have plotted in Fig. El the level scheme for the H 2 and 
H 2 systems. We consider first the situation with only one ls-like orbital per atom. 
For H 2 we have („) = 6 two particle states For that purpose, we start with the 
parameterized Hamiltonian (J60|) . where subscripts "i" and "j" label now the two atomic 
sites and hence U\ = U 2 = U, K X2 = K, J 12 = J, V\ 2 = V, and e\ = e 2 = e a . The lowest 
eigenstate for H 2 is 

E G = A 5 = 2e a + 1 -{U + K) + J - l - [(U - Kf + 16(t + V) 2 ] 1/2 , (61) 
and the corresponding singlet ground state in the Fock space has the form 

\G>= 1 x 

y/2D{D -U + K) 

{ 4 ^/2^ (°1T°L ~ a U a 5t) ~ ~ ^ + ~ ( at it a 2| + a u at 2T)} 1° >> ( 62 ) 
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where 

D= [(U -K) 2 + lQ{t + V) 2 ] 1/2 . 

The lowest spin-singlet eigenstate has an admixture of symmetric ionic state ^(a^a^-l- 
aj.aL). Therefore, to see the difference with either the Hartree-Fock or Heitler-London 
approach to H 2 we construct the two-particle wave function for the ground state 
according to the prescription 

$ (ri,r 2 ) = < 0|*( ri )#(r 2 )|G > . (63) 
Taking *(r) = ELi Ei= T M r )x*{*), we obtain that 



2(t + V) ^ , , 1 D-U + K^ , . 
$ n,r 2 = ^ e (ri,r 2 ) — \/ $;(ri,r 2 ), 64 

where the covalent part is 

$ c (ri,r 2 )= (65) 

[w 1 (r 1 )w 2 (r 2 ) + twi(r 2 )tw 2 (ri)] [xi(*i)Xl( r 2) ~ Xli^Xlfa)] , 
whereas the ionic part takes the form 

$i(ri,r 2 )= (66) 

[wi(ri)wi(r 2 ) + w 2 (r 1 )w 2 (r 2 )] [Xt( r i)Xl( r 2) - X|(ri)XT( r 2)] • 

The ratio of the coefficients before $ c (ri,r 2 ) and $j(ri,r 2 ) can be regarded as the 
many-body covalency 7^. This value should be distinguished from the single-particle 
covalency 7 appearing in the definition of the orthogonalized atomic orbital u>j(r): 

<r)=/3[$,,(r)-7$^r)], (67) 

with j 7^ i. The two quantities are drawn in Fig. El The many-body covalency 7 m b 
represents a true degree of multiparticle configurational mixing. 

In Table El we list the energies and the values of the microscopic parameters for 
H 2 system with optimized orbitals, whereas in Table 0] the same is provided for 
molecular ion. One should notice a drastic difference for the so-called correlated hopping 
matrix element V in the two cases. The same holds true for the direct exchange integral 
J (ferromagnetic). This exchange integral is always decisively smaller than that for the 
antiferromagnetic kinetic exchange, J^ ex = 4(t + V) 2 / (U — K). The virtual interatomic 
hopping processes leading to the strong kinetic exchange are the source of the singlet 
nature of the H 2 ground state. The H 2 ground state is unstable with respect to the 
dissociation into H 2 and e~, contrary to the H~ case. However, the energetics of such 
state is important when calculating e.g. the metallization of molecular hydrogen or 
determining the tunnelling conductivity through H 2 molecule, as shown schematically 
in Fig. El This last Figure illustrates the method of determining the energetics of excited 
states of H 2 by measuring e.g. the tunnelling conductivity (i.e. via H 2 intermediate 
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2 4 6 8 10 

interatomic distance R/^ 



Figure 3. The single-particle (7) and many-body (^ m b) covalency factors for the H2 
wave functions. For details see main text. 

state, see below). 

For the sake of completeness, we have provided in Fig. 01 the energy levels vs. R 
for H21H2 , and ions and have labeled both the bonding (B) and antibonding (AB) 
level positions. We see that the distance between H 2 and levels is much smaller 
then the distance between H 2 and H 2 states. This will lead to the asymmetry in the 
tunneling conductivity when reversing the bias voltage sign, as discussed in detail in 
Section 5. 

3.4- Hydrogen clusters, 

As next application we consider hydrogen-cluster systems, with N < 6 atoms. We 
take the atomic-like Is orbitals {$j(r)} of an adjustable size a = a -1 , composing the 
orthogonalized atomic (Wannier) functions {wi(r)} i= i ^. The cluster of TV atoms with 
N electrons contains ( 2 ^) states and the second-quantized Hamiltonian is of the form 
with 3- and 4- site terms added. The 3- and 4- site interaction terms are difficult 
to calculate in the Slater basis (see below). Therefore, we have made an ansatz [7] 
namely, we impose the condition on the trial Wannier function that the 3- and 4-site 
matrix elements vanish. This allows for an explicit expression of the 3- and 4- 
site matrix elements V( jkl in the atomic representation via the corresponding one- and 
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interatomic distance, R (a B ) 



Figure 4. Position of H% bonding(B) and antibonding (AB) energy levels with respect 
to those of H 2 and Ho, drawn as a function of intermolecular distance. 



Table 3. Ground-state energy and microscopic parameters (in Ry) for H 2 molecule. 
The last column represents the kinetic exchange integral characterizing intersite 
antifcrromagnetic exchange 



R/a 


E G /N 


e« 


t 


U 


K 


V [mRy] 


J [mRy] 




1.0 


-1.0937 


-1.6555 


-1.1719 


1.8582 


1.1334 


-13.5502 


26.2545 


7755.52 


1.5 


-1.1472 


-1.7528 


-0.6784 


1.6265 


0.9331 


-11.6875 


21.2529 


2747.41 


2.0 


-1.1177 


-1.722 


-0.4274 


1.4747 


0.7925 


-11.5774 


16.9218 


1130.19 


2.5 


-1.0787 


-1.6598 


-0.2833 


1.3769 


0.6887 


-12.0544 


13.1498 


507.209 


3.0 


-1.0469 


-1.5947 


-0.1932 


1.3171 


0.6077 


-12.594 


9.8153 


238.939 


3.5 


-1.0254 


-1.5347 


-0.1333 


1.2835 


0.5414 


-12.8122 


6.9224 


115.143 


4.0 


-1.0127 


-1.4816 


-0.0919 


1.2663 


0.4854 


-12.441 


4.5736 


55.8193 


4.5 


-1.006 


-1.4355 


-0.0629 


1.2579 


0.4377 


-11.4414 


2.8367 


26.9722 


5.0 


-1.0028 


-1.3957 


-0.0426 


1.2539 


0.3970 


-9.9894 


1.6652 


12.9352 


5.5 


-1.0012 


-1.3616 


-0.0286 


1.2519 


0.3623 


-8.3378 


0.9334 


6.1455 


6.0 


-1.0005 


-1.3324 


-0.01905 


1.251 


0.3327 


-6.7029 


0.5033 


2.8902 


6.5 


-1.00024 


-1.3073 


-0.0126 


1.2505 


0.3075 


-5.2242 


0.2626 


1.3452 


7.0 


-1.0001 


-1.2855 


-0.0083 


1.2503 


0.2856 


-3.9685 


0.1333 


0.6197 


7.5 


-1.00004 


-1.2666 


-0.0054 


1.2501 


0.2666 


-2.9509 


0.066 


0.2826 


8.0 


-1.00002 


-1.25 


-0.0035 


1.25006 


0.25 


-2.1551 


0.032 


0.1277 


8.5 


-1.00001 


-1.2353 


-0.0023 


1.25003 


0.2353 


-1.5501 


0.01523 


0.0572 


9.0 


-1. 


-1.2222 


-0.0015 


1.25001 


0.2222 


-1.1005 


0.0071 


0.0254 


9.5 


-1. 


-1.2105 


-0.0009 


1.25001 


0.2105 


-0.7725 


0.0033 


0.0112 


10.0 


-1. 


-1.2 


-0.0006 


1.25 


0.2 


-0.5371 


0.0015 


0.0049 
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Table 4. Same as in Table |31 for H 2 ion 



R/ 


a 


E G /N 




e a 




t 




u 




K 




V 




J 


A(t+vy 

U-K 


1. 


.0 


-0 


.4591 


-1 


.6607 


-0. 


.5869 


1 


.1414 





.7360 


-0 


.0105 


0. 


.0163 


3.5220 


1. 


.5 


-0 


.7659 


-1 


.6647 


-0. 


,4285 


1 


.1279 





.6983 


-0 


.0085 


0. 


.0161 


1.7782 


2. 


.0 


-0 


.8813 


-1 


.6259 


-0. 


,3083 


1 


.0979 





.6474 


-0 


.0078 


0. 


.0150 


0.8871 


2. 


.5 


-0 


.9264 


-1 


.5737 


-0. 


,2221 


1 


.0692 





.5961 


-0 


.0079 


0. 


.0133 


0.4476 


3. 


.0 


-0 


.9423 


-1 


.5204 


-0. 


,1603 


1 


.0466 





.5476 


-0 


.0086 


0. 


.0113 


0.2286 


3. 


.5 


-0 


.9460 


-1 


.4704 


-0. 


.1154 


1 


.0305 





.5025 


-0 


.0093 


0. 


.0091 


0.1179 


4. 


.0 


-0 


.9450 


-1 


.4252 


-0. 


.0826 


1 


.0196 





.4608 


-0 


.0099 


0. 


.0071 


0.0612 


4. 


.5 


-0 


.9426 


-1 


.3848 


-0. 


.0585 


1 


.0126 





.4226 


-0 


.0101 


0. 


.0052 


0.0319 


5. 


.0 


-0 


.9402 


-1 


.3491 


-0. 


.0410 


1 


.0080 





.3881 


-0 


.0099 


0. 


,0037 


0.0167 


5. 


.5 


-0 


.9384 


-1 


.3176 


-0. 


.0284 


1 


.0051 





.3573 


-0 


.0093 


0. 


.0025 


0.0088 


6. 


.0 


-0 


.9373 


-1 


.2901 


-0. 


.0194 


1 


.0032 





.3300 


-0 


.0085 


0. 


.0017 


0.0046 


6. 


.5 


-0 


.9365 


-1 


.2621 


-0. 


.0130 





.9905 





.3058 


-0 


.0075 


0. 


.0011 


0.0025 


7. 


.0 


-0 


.9363 


-1 


.2402 


-0. 


0086 





.9876 





.2847 


-0 


.0065 


0. 


,0007 


0.0013 


7. 


.5 


-0 


.9365 


-1 


.2211 


-0. 


.0056 





.9856 





.2662 


-0 


.0055 


0. 


.0004 


0.0007 


8. 


.0 


-0 


.9367 


-1 


.2044 


-0. 


.0036 





.9844 





.2498 


-0 


.0046 


0. 


.0003 


0.0004 


8. 


.5 


-0 


.9372 


-1 


.1897 


-0. 


,0022 





.9839 





.2352 


-0 


.0037 


0. 


0002 


0.0002 


9. 
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.9376 


-1 


.1768 
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.2222 


-0 


.0030 
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0.00010 


9. 
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-0 
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-1 
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-0. 


,0008 
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-0 


.0024 
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10. 


.0 


-0 


.9384 


-1 


.1549 


-0. 
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.9848 





.2000 


-0 


.0018 


0.00003 


0.00003 



two- site elements. In Figs. Eland El we present the results for the ground- and excited- 
states energies for the square and face-centered-square (fcs) configurations, N = 4 and 
5, respectively. The states are grouped into manifolds, which are characterized by 
the number of double occupancies of a single state u>i(r), appearing in the system. 
The horizontal lines mark the ground state, states with one and two double electron 
occupancies in the atomic limit (i.e. for large interatomic distance). The manifolds thus 
correspond to the Hubbard subbands introduced for strongly correlated solids [2E]- As 
far as we are aware of, our results are the first manifestation of the energy manifold 
evolution into well separated subbands with the interatomic distance increase. The 
first two subbands correspond to HOMO and LUMO levels determined in quantum- 
chemical calculation [20] • In. Fig. Q and El we represent the renormalized Wannier 
function profiles for the face centered square configuration of N = 5 atoms, for the 
central and the corner positions, respectively. Note the small negative values on the 
nearest-neighbor sites to assure the orthogonality of the functions centered on different 
sites. Obviously, the atom in the center of the square is inequivalent to the remaining 
four corner atoms, as can be see explicitly in Fig. |H1 where the density profile n(r) in 
the cluster plane, according to formula (@7j), has been drawn. The electron density near 
the central atom is decisively smaller, a clear sign of electron-correlation effects induced 
by the repulsive Coulomb interaction. This assessment is corroborated in Fig. UHl 
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Figure 5. Ground- and excited states energies for the H4 square configuration as a 
function of the interatomic distance. 

where the difference to the Hartree-Fock part of the density profile has been presented. 
These densities should be possible to be determined with the help of scanning tunneling 
microscopy (STM). 

3.5. Energetically stable H4 clusters and fermionic nanoladders 

The Hn clusters of regular-polygon shape are not stable energetically as the ground 
state energies in Figs. El and El at the minimum (per atom) is above that for the H2 
case. This is because a stable, say H4 cluster will reflect strong molecular bonding of 
each pair of H 2 molecules to saturate the covalent character of the bonding. To prove 
that this is indeed the case we have considered an example of a rectangular cluster 
differentiating between the lateral distance (say along the bond of the length a) and the 
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Figure 6. Ground- and excited states energies for the face-centered-square (fcs) 
configuration as a function of the interatomic distance. 

horizontal intermolecular distance b. We consider both the planar and the twisted by 
90° configurations, as drawn in Fig. [TT] The Gaussian STO-3G basis of an adjustable 
size a~ l has been used in this case, so the 3- and 4- site interaction terms in the atomic 
basis are included explicitly. Additionally, the Hamiltonian in the second quantization, 
only containing the principal one- and two- site interactions, is taken in the form [S] 

H = e a ff Hia + % a L°W + U n ^ ni i + \ K ij Sn i$ n j, ( 68 ) 

where 5ui = 1 — n i; rii = J2 a n io . and 
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Figure 7. The renormalized Wannier function profile for for the central atom in the 
fcs structure. 



with the last term in the parenthesis representing the Coulomb repulsion (in atomic 
units) between the protons placed at the distance Ry, together with the electron-electron 
intersite repulsion (such a redefinition is necessary to achieve the proper atomic limit, 
when the nearest-neighbor distance R — ► oo. Also, we have to take into account three 
different hopping integrals: between the nearest neighbors t\ = ti 2 = t 34: , between the 
second neighbors t 2 = ti3 = £24, and the third neighbors £3 = tu = £23- Likewise, we 
have three intersite Coulomb interactions {i^ n }n=i,2,3- 

The basic ground-state characteristics for the two spatial arrangements is provided 
in Tables El and El 

The results are listed as a function of intermolecular distance and include: the bond 
length a m i n , the inverse Gaussian-function size a, the energy of H± and the energy of 
two H 2 molecules (per atom). The principal feature of these results is that the system 
is energetically stable against the dissociation into two H 2 molecules (cf. the last row 
in the two situations). The global minima are: 

i) for the planar geometry: b min = 4.32a , a min = 1.4303a , ot min = 1.1937oq , and 
E ™in/ N = -1.149061%; and 

ii) for the twisted 90° geometry: b min = 4.13ao, a m in — 1.4318ao, «mm = 1.11927ao , 
and E™ in /N = -1.150396%. 
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Figure 8. The rcnormalized Wannier function profile for for the corner atom in the 
fcs structure. 



Table 5. The optimal bond length a m i n , inverse orbital size a m ; n , and the ground- 
state energy Eg per atom for the planar H4 cluster (i). The corresponding energy of 
the molecular dimer Eg (2H2) is also provided. 



b/a 


<3min/ do 




E G /N 


E G (2E 2 )/N 


1.7 


1.3627 


1.2231 


-0.928424 


-0.922411 


2.0 


1.3291 


1.2395 


-1.019152 


-1.016365 


2.5 


1.3518 


1.2304 


-1.098551 


-1.097314 


3.0 


1.3829 


1.2157 


-1.131770 


-1.131167 


3.5 


1.4075 


1.2041 


-1.144613 


-1.144317 


4.0 


1.4238 


1.1969 


-1.148598 


-1.148454 


5.0 


1.4373 


1.1911 


-1.148093 


-1.148056 


6.0 


1.4390 


1.1908 


-1.145975 


-1.145964 


8.0 


1.4375 


1.1924 


-1.143651 


-1.143649 


10.0 


1.4366 


1.1929 


-1.142756 


-1.142755 


20.0 


1.4357 


1.1940 


-1.141908 


-1.141908 


00 


1.4356 


1.1943 


-1.141783 
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Figure 9. Electron density profile for the H§ face-centered-square (fcs) structure. 



So, the geometry (ii) is the most stable in vacuum and should constitute a building 
block for i?2 molecular crystal. Note we have not included the zero point motion of the 
nuclei j2]. The ratio b/a is roughly four, so substantial molecular identity of H 2 pairs 
survives. 

As in all cases before, the knowledge of microscopic parameters is crucial for the 
larger-cluster or solid-state configurations ( Hamiltonian ()68|) provides whole dynamics 
within the subspace with one orbital per atom). Therefore, in and El we list them for 
the two configurations considered, as a function of intermolecular distance b ( the last 
row represents the corresponding values for H 2 in the ground state). 

The stability of the H± cluster raises a very interesting question of stability of H 2 n 
ladders, which can be regarded as the simplest fermionic ladders with the frustration of 
electron spins in the twisted configuration. 
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Figure 10. Difference between the calculated density of states and the Hartree-Fock 
results for the H$ fcs structure. 




Figure 11. Schematic representation of the H4 cluster geometries: (i) parallel and 
(ii) perpendicular orientation of H2 molecules. Geometrical parameters of the cluster 
are: the bond length a and the intermolecular distance b. The numbering order of the 
lattice sites j is also provided. 
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Table 6. The optimal bond length a m i n , inverse orbital size a m ; n , and the ground- 
state energy Eq per atom for the cluster geometry (ii). The corresponding energy of 
the molecular dimer Eg(2H2) is also provided. 



b/a 


a min/ a 




E G /N 


E G (m 2 )/N 


1.6 


1.5796 


1.1568 


-0.928235 


-0 


.923390 


2.0 


1.3759 


1.2206 


-1.038891 


-1 


.037968 


2.5 


1.3725 


1.2193 


-1.108116 


-1 


.107803 


3.0 


1.3947 


1.2091 


-1.136931 


-1 


.136809 


3.5 


1.4153 


1.2000 


-1.147566 


-1 


.147519 


4.0 


1.4292 


1.1937 


-1.150310 


-1 


.150293 


5.0 


1.4397 


1.1894 


-1.148674 


-1 


.148672 


6.0 


1.4400 


1.1897 


-1.146200 


-1 


.146200 


8.0 


1.4377 


1.1926 


-1.143705 


-1 


.143705 


10.0 


1.4367 


1.1929 


-1.142774 


-1 


.142774 


20.0 


1.4357 


1.1940 


-1.141908 


-1 


.141908 


oo 


1.4356 


1.1943 


-1.141783 







Table 7. Microscopic parameters (in Ry) for the H4 cluster configuration (1) calculated 
in the Gaussian STO-3G basis. Corresponding values of the inverse orbital size a m i n 
and the bond length a m j n are provided in Table 



b/a 


ffiff 






t 2 




*3 


u 


K x 


K 2 


K 3 


1.7 


-0. 


.2354 


-0. 


,8610 


-0, 


,6622 


0.0822 


1. 


,811 


1 


.027 


0. 


,947 


0.728 


2.0 


-0. 


,3088 


-0, 


,8791 


-0, 


,5137 


0.0617 


1. 


,802 


1 


.032 


0, 


,872 


0.692 


2.5 


-0. 


,4233 


-0. 


,8390 


-0. 


,3268 


0.0352 


1. 


,748 


1 


.007 


0. 


,750 


0.619 


3.0 


-0. 


,4925 


-0. 


,7983 


-0. 


,2067 


0.0220 


1. 


,702 





.984 


0. 


,649 


0.553 


3.5 


-0. 


,5319 


-0, 


,7685 


-0, 


,1290 


0.0150 


1, 


,671 





.967 


0, 


,567 


0.497 


4.0 


-0. 


5533 


-0, 


,7492 


-0, 


,0785 


0.0103 


1, 


653 





.957 


0, 


,500 


0.449 


5.0 


-0. 


,5689 


-0, 


,7316 


-0, 


,0267 


0.0037 


1, 


639 





.949 


0, 


,401 


0.373 


6.0 


-0. 


,5707 


-0. 


,7275 


-0. 


.0084 


0.0007 


1. 


,638 





.948 


0. 


,334 


0.318 


8.0 


-0. 


,5692 


-0. 


,7269 


-0. 


,0006 


-0.0000 


1. 


,640 





.949 


0. 


,250 


0.243 


10.0 


-0. 


,5684 


-0, 


,7268 


-0, 


,0000 


-0.0000 


1, 


,641 





.950 


0, 


,200 


0.197 


20.0 


-0. 


5673 


-0, 


,7272 


-0, 


,0000 


-0.0000 


1, 


,642 





.951 


0, 


,100 


0.100 
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Figure 12. Ground-state energy per atom for the planar ladder of TV =12 atoms (a) 
with a fixed intermolecular distance b (specified in the units of Bohr radius ao), and (b) 
the corresponding values of the optimal inverse orbital size a m in- Periodic boundary 
conditions are used. The results for the atomic chain (dashed lines) with a lattice 
parameter a < 6 C rit are shown for the comparison. 



For that purpose we have considered planar (i) and 90° twisted (ii) ladder composed 
of 8 4- 12 atoms. In Fig. ^] we display the ground state energy (per atom) for 
the nanoladder containing N = 10 atoms with periodic boundary conditions. The 
insets provide the values of the optimal bond length a min (left) and the inverse atomic 
(Gaussian) orbital size a min (right). The Gaussian orbitals (STO-3G basis) have been 
used to define the single-particle basis. The horizontal dashed line marks H2 ground 
state energy. One should also note that due to the closed-shell molecular-crystal 
configuration (b/a ~ 4) of the ground state, the data almost do not depend on the 
system size (e.g. analogical results for iV = 8 fit onto those shown in Fig. H2]up to the 
pixel size). The characteristics of the global minima are: 

i) for the planar geometry: b min = 4.00a , a m i n = 1.422a , ct m i n = 1.189cIq , and 
E™ in /N = -1.1626%; and 
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Figure 13. Phase diagram of the planar fermionic ladder. The optimal value of the 
bond length a m j n is drawn as a dashed line. Dotted line marks the effective quarter- 
filling obtained from the noninteracting band model. 

ii) for the twisted 90° geometry: b min = 3.67a , a min = 1.426a , ce min = 1.178oq , and 
E™ in /N = -1.1680%. 

These values lead to the binding energies of the molecules in the nanoladders: 
AEq /N = 20.8mRy and AEq 1 ^ /N = 2Q.2mRy, respectively, i.e. about three times 
larger than for the H± clusters. 

One can also draw some conclusions about the nature of electronic states in those 
nanoladders. Namely, by calculating the so called charge stiffness (Drude weight) and 
the charge gap (cf. Ref. [3]) one can draw a phase diagram shown in Fig. 13 for a planar 
nanoladder on the plane a — b. For either a or b large we expect the Mott insulating state 
of spins S = 1/2. However, it is interesting to note that for bond lengths in the range 
2-^3 and an appropriate intermolecular distance (b ~ a) a quasimetalic state is possible 
and is followed by a band insulator for small b < a. The metal phase sandwiching the 
two, band and Mott insulating phases appears well beyond the optimal configuration 
for the nanoladder (see the dashed line representing a = a m i n as a function of b). So it 
can appear only when nanoladder is artificially made on a supporting substrate. The 
dotted line marks the model of qurted-filled (QF) band system discussed in Ref. [5]. 
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4. Monoatomic nanowires: from nanometal to Mott insulator 

4-1. Ground state properties 

The EDABI method was applied first to nanochains starting with Slater basis of atomic 
states CHI EEZj In this Section we present the new results starting from STO-3G 
"atomic" basis of adjustable size. The selection of this basis makes the inclusion of 3- 
and 4-site interactions possible. We limit here to the situation with one electron per 
atom. Strictly speaking, we study a nanochain of hydrogen atoms, which can also model 
the behavior of single valence electron nanowires. 

We start again with Hamiltonian (ffiHJ) describing the so-called extended Hubbard 
model for a system with one orbital per atom (here the orbitals are taken as Gaussians of 
an adjustable size, out of which we compose the Wannier functions). The diagonalization 
in the Fock space is performed with the help of Lanczos method described in detail 
elsewhere j3 E|. In Fig. H3I we plot the ground state energy of chains containing 
N = 6 — 10 hydrogen atoms, as a function of interatomic distance R (in units of 
do). The Wannier functions are calculated in the tight-binding approximation, with six 
attractive atomic Coulomb wells representing the periodic potential. In the inset we 
display the inverse size a (in units of inverse ao) of the Gaussian functions (note that 
the actual size of atoms is reduces in the correlates state). The continuous Eg lines INS 
and M represent respectively to the Mott insulating state represented by E G NS = e e Jl ', 
and to the ideal metallic state, for which 



For the sake of completeness, we have also included the Hartree-Fock (HF) result for 
Eg (the dotted line) calculated for the antiferromagnetic (Slater) state; it provides a 
better estimate of Eg than either M or INS curves. For discussing HF solution we take 
the Hamiltonian (|68|) in the form 




(70) 



where the correlation is taken for the one-dimesional ideal gas 




(71) 






i<j 




(72) 



The combined exact diagonalization- ab initio approach 



34 



Table 8. Microscopic parameters (in Ry) for the H4 cluster configuration (ii) 
calculated in the Gaussian STO-3G basis. Corresponding values of the inverse orbital 
size a m i n and the bond length a min are provided in Table El 



b/a 


ffiff 


t x 


*2 


h 


u 


K x 


K 2 


K 3 


1.6 


-0.5352 


-0.5423 


-0.2986 


-0.2986 


1.650 


0.911 


0.816 


0.816 


2.0 


-0.4511 


-0.7496 


-0.2230 


-0.2230 


1.732 


0.991 


0.772 


0.772 


2.5 


-0.4863 


-0.7820 


-0.1451 


-0.1451 


1.711 


0.986 


0.680 


0.680 


3.0 


-0.5221 


-0.7721 


-0.0923 


-0.0923 


1.683 


0.972 


0.598 


0.598 


3.5 


-0.5467 


-0.7558 


-0.0571 


-0.0571 


1.661 


0.961 


0.530 


0.530 


4.0 


-0.5613 


-0.7423 


-0.0342 


-0.0342 


1.646 


0.953 


0.473 


0.473 


5.0 


-0.5715 


-0.7293 


-0.0115 


-0.0115 


1.637 


0.947 


0.387 


0.387 


6.0 


-0.5719 


-0.7264 


-0.0038 


-0.0038 


1.636 


0.947 


0.326 


0.326 


8.0 


-0.5692 


-0.7269 


-0.0003 


-0.0003 


1.640 


0.949 


0.247 


0.247 


10.0 


-0.5685 


-0.7268 


-0.0000 


-0.0000 


1.641 


0.950 


0.198 


0.198 


20.0 


-0.5673 


-0.7272 


-0.0000 


-0.0000 


1.642 


0.951 


0.100 


0.100 


00 


-0.5671 


-0.7273 
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Figure 14. Ground state energy per atom vs. R for the linear chain with N — 
6 -j- 10 atoms with periodic boundary conditions. The STO-3G Gaussian basis for 
representation of atomic orbitals forming the Wannier function has been used. The 
inset provides a universal behavior of the inverse size a of the orbitals. For details see 
main text. 
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Figure 15. Parity effect on spin ordering: spin-spin correlations for nanochains of 
N = 10 (a) and N — 11 (b) atoms. The values of the interatomic distance R are 
specified in the atomic units (do = 0.529 A). 

The HF solution involves calculations of the ground-state energy with a simultaneous 
self-consistent determination of the sublattice magnetization m = (n^ — n^) and of 
the hopping correlation function (a^a^), which will not be discussed in detail here |S] 
(obviously, (Srii) =0). 

The most spectacular are the spin-spin correlations in the collective spin-singlet 
(YliLi = 0) state- in Fig- IT3k we display the corresponding spin-spin correlation 
function (Sj ■ Sj) as a function of the distance — for N = 10 and several values of R. 
We observe Neel-like state correlations in this spin singlet state as it oscillates through 
the whole system length. For N = 11 these quasi oscillatory correlations are strongly 
suppressed due to the spin-frustration effects, as shown explicitly in Fig. ITHb . Changing 
the lattice constant R does not alter the picture qualitatively. This long-range feature of 
spin-spin correlations will have profound consequences on the single-particle spectrum, 
as we discuss below. Here the boundary conditions discussed below (cf. Sec. 4.5) have 
been considered, although the difference between the and b survives even for the 

periodic boundary conditions. Additional features of the N dependence of < Sj ■ Sj > 
have been discussed in j2H 

4-2. Momentum distribution in the ground state: Fermi or Luttinger nano liquid 
localization? 

We now introduce the particle quasimomentum distribution in the ground state for 
the situation with one electron per atom; this is displayed in Fig. EH for the number of 
atoms iV = 6 -f- 14. We observe a very universal character of the curves provided the 
periodic boundary conditions are taken for the chains of N = 4n + 2 atoms, whereas 
the antiperiodic boundary conditions are taken for N = 4n atoms, with n being a nat- 
ural number. It is very tempting to regard the distributions for R/a < 3.5 as modified 
Fermi-Dirac function characterizing the Fermi liquid even for such short chains for which 
we have discrete momentum states. However, for R/a Q — > 4 the distributions becomes 



The combined exact diagonalization- ab initio approach 



36 



< 

OL- 
ID 
O 

o 
o 
z 
o 

DC 
h- 

o 

LU 
_l 
LU 







R/a =1,5^ +-o **<o i o* 



Q>« >KOi 



0.8 
0.6 
0.4 
0.2 Jf >» 



R/a =2.5* +• 

R/a =3.5^+?? ; „ r 
& R/a =4.5 



,o * * o, 



N=6 
N=8 
x N=10 
\ N=12 
•+ * N=14 



x 

G 



4 o>#^u 

0>«>KO«+>^ 



'O 

o * »<o -i 



-0.5 0.5 

MOMENTUM, kR/Jt 



Figure 16. Momentum distribution nko- for the number of atoms N — 6 -j- 14 as it 
evolves with the increasing lattice constanti?. 
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Figure 17. Statistical distribution for electrons in a chain of N = 6 -j- 14 atoms 
with periodic boundary conditions. The interatomic distance i? is specified in units of 
Bohr radius a$. The continuous line represents the parabolic interpolation, which is of 
the same type for both k > &f and k < fop. 





Figure 19. Luttinger-liquid scaling for a half-filled one-dimensional chain of N = 
6 -j- 14 atoms. Interatomic distance R is specified in the units of a . 
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Table 9. Optimized inverse orbital size, microscopic parameters and the ground-state 
energy for N — 10 atoms calculated in Slater-type basis, as a function of interatomic 
distance. Intersite Coulomb repulsion K\ is included on the mean-field level in e e Jf , 
Hubbard U term is treated exactly. Single-particle potential contains six Coulomb 
wells. 
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Figure 20. Spectral functions A(k, u) for N = 10 (a) and N — 12 atoms, with periodic 
boundary conditions. A clear dispersion of the states is observed for 7? = 1.5ao, which 
transform into the atomic peaks with the increasing R, through an incoherent regime 
for R/a n fa 3 
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Figure 21. Density of single-particle states for N = 10 atoms. Note a pronounced 
incoherent spectrum for R in the range 3 -j- 5ao 




Figure 22. Spectral-density-peak positions for the nanochain of N = 10 (left panel) 
and N = 11 (right panel) atoms with generalized boundary conditions. The Hartree- 
Fock (solid line) and noninteracting system (dashed line) dispersion relations are shown 
for comparison. 
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continuous, i.e. without an apparent Fermi ridge at Fermi momentum kp = tt/(2R), in 
rough agreement with the preliminary results for the case with the Slater-type orbitals 
[T3~] . So, a crossover with the increasing interatomic distance from the Fermi liquid-like 
behavior to the chain of atomic-like states is clearly seen. It should be noted that the 
Luttinger liquid-like scaling fitting to the points displayed in Fig. EI] is equally convinc- 
ing, at least for this half-filled band configuration, as discussed next. 

4-3. Single-particle spectral density function evolution 

The state of a nanochain viewed through the statistical distribution of momentum 
states can be discussed as follows. The ideal Mott localized state, as is seen in the 
bulk systems, is not possible for finite system. This is simply because the length of 
the system is finite and therefore, the probability of the electron tunneling from an 
atomic state on one chain end to the other end is nonzero. The quantitative question 
is, whether it is of the order e~ NaR or larger. Due to the fact that the overlap integral 
is ~ ^(aR) 2 e~ aR , the probability is enhanced by a nonzero hopping amplitude between 
more distant neighbors. Having in mind this circumstance, it is at least imaginable, that 
the distribution function displayed in Fig. EI] is Fermi-Dirac like for aR ~ 1. One can 
thus try to formalize this observation further. Namely, in Fig. El we plot for a half 
filled band system containing up to iV = 14 electrons. The continuous lines represent 
the parabolic parametrization: 

ft-ka — -sgn(k — kp) [a(k — kp) 2 + f3\k — kp\ + 7] , (73) 
Table 10. Same as in Tabled] calculated in Gaussian-type STO-3G basis. 
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for both k < kp and k > kp (note that the Fermi point is not occupied for N even). This 
parametrization, in fact implies a discontinuous distribution, allows us to interpret the 
distribution discontinuity Arik F = rik F -o — n k F +o in terms effective mass enhancement 
at the Fermi point. 

m F /m BAND = (ArikpY 1 = Z' 1 , (74) 

where Z~ l is the usual Fermi-liquid enhancement factor. The corresponding 
enhancement factor determined in that manner is plotted in Fig. (squares), whereas 
the solid line represents the finite-size scaling with interatomic distance, mp / m ban d = 
A\R — Rq^ 1 i with A ~ 10.2, the localization thereshold is Redo = 3.92, and the critical 
exponent has the approximate value 7 = 4/3. This behavior emulates a quantum critical 
behavior and should not be taken literally. Nonetheless, the value of Rc distinguishes 
qualitatively between the nanometallic state R < Rc (for which the distribution (|75|l has 
a jump at hp) and the Mott insulating (semiconducting) state (R > Rc, the distribution 
function is a continous function of k). 

However, the situation is not that simple. One can ask the question whether this 
short one-dimensional system is not rather resembling the Luttinger liquid-like behavior. 
This task has been undertaken seriously and in Fig. EI] we plot the fitting of the data 
displayed in Fig. ITBlto the corresponding dependence for the Tomanaga-Luttinger liquid, 
with the logarithmic correction included, namely [27] 

In \nk F — rik(j\ = — © In z + 6 In In z + c + 9(1/ In z), (75) 

where z = ir/\kp — ir/R\,Q,b, and c are constants. The parameters depend on the 
distance R, as shown elsewhere j^Z] • This fitting provides also the localization thereshold 
Rc — 2.60ob, which is reduced by 50% from the value for the Fermi liquid interpretation 
of Uka- One should also note that the present interpretation of the exact solution 
displayed in Fig. HB] does not admit a discontinuity at k = kp. However, it is quite 
amazing that both the Fermi- and Luttinger-liquid interpretations can provide satisfying 
interpretation of n^ a to an equal degree. This means that there must be an underlying 
universal behavior of a new type, incorporating Fermi- and Luttinger-liquid concepts, 
at least semiquantitatively. Nonetheless, for the not too large R the Fermi-Dirac-like 
distribution fits better the k dependence. The situation for odd iV requires an explicit 
discussion of boundary conditions and is provided below. 

To characterize the spectrum of single-particle excitations we use the definition of 
the spectral-density function 

A(k,a;) = El^k^f 5 [u> - (E^ - <)] , (76) 

n 

where upper (lower) sign corresponds to the energies with u > \i (u < /i), respectively, 
|3>„ ) is the n-th eigenstate of the system containing N particles, is the corresponding 
eigenenergy, and < .. > is the matrix element are calculated within the Lanczos 
algorithm (221 • l n Fig. HHl we present the panel of A(k,uj) for three distances R. For 
small R, a clear two-peak structure appears at the Fermi momenta kp = ±n/(2R) for 
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N = 12; an artificial broadening of the peaks appears because we use the approximate 
expression for the S(x) function: 5(x) ~ (l/7r)e/(x 2 + e 2 ), with e = 10~ 2 Ry. The 
splitting is caused by the antiferromagnetic correlations depicted in Fig. ITHk . In the 
range R/oq < 2 the quasiparticles are well defined, but an incoherent part grows with the 
increasing R. For R/ao ~ 3-^4 the Hubbard subbands are formed and evolve continuously 
into atomic levels located respectively at u = e a and lj = e a + U for R — ► oo. Those two 
limiting peak positions correspond to the ground H° and excited H~ states, respectively. 
Combining the last results with the corresponding discussion for the Hn clusters one sees 
that the Hubbard subband structure represents a universal feature of nanoscopic systems. 
In the limit of larger interatomic separation this structure is clearly distinguishable from 
the discrete level structure coming from geometric quantization of this confined system. 

To demonstrate the importance of incoherent part of the spectrum we have 
calculated the density of states Af(u) = ^ k A(k, a;), displayed in Fig. ED for N = 10 
atoms, for the interatomic distances specified. In the regime R = 3 4- 5a well resolved 
quasiparticle peaks coalesce into a complicated random-like (incoherent) spectrum, out 
of which clean atomic peaks emerge for larger R. Note that the intermediate regime 
corresponds to the situation for which the bare bandwidth W = 4|t| of the single particle 
states fulfils the condition W ~ U, i.e. the electronic system switches from the weak- to 
the rather strong-correlations regime, corresponding to the delocalization-localization 
crossover of the Mott-Hubbard type taking place. For one electron per atom the lower 
energy manifold is filled with electrons, whereas the upper is empty. The presence of 
the incoherent spectrum for R « 3 -f- 5ao seems to represent also a universal feature, as 
it appears also for larger values of N. 

4-4- Slater vs. Gaussian trial basis 

We now compare the results obtained within EDABI when either the adjustable Slater 
or the STO-3G functions are used. Probably, most interesting is to compare the results 
for the Hubbard model with intersite Coulomb interactions taken in the Hartree-Fock 
approximation (only then the results for the Slater ls-like basis are energetically stable. 
In Table El we list the results for the Slater basis containing both Eq and microscopic 
parameters for a ring of iV = 10 atoms. The same for the case of STO-3G basis is shown 
in Table One sees that the Gaussian basis provides a lower value of Eq. This is 
because in the Slater case we have neglected 3- and 4- site interactions in the atomic 
basis that represents a rather crude approximation. On the contrary, the results for 
Eq obtained for the STO-3G basis for N = 10 atoms are very close (within 10~ 2 Ry) 
to those obtained for an exact solution for an infinite chain [3U| and the values of the 
microscopic parameters within 10~ 3 [2E]- Because of the good accuracy of the results for 
the Gaussian basis we have also provided in Table ITT1 the results for the chain of iV = 10 
atoms when the long-range part of the Coulomb interactions ~ are included. These 
results provide the range of variation of the microscopic parameters when compared to 
either H 2 (cf. Table |2J) or U\ value for He atom (cf. Table |3J). 
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4-5. The role of boundary- conditions: parity effects 

It is believed that system properties in the thermodynamic limit (N — > 00) are the same 
regardless of the boundary conditions used. This claim has been tested on many model 
systems. However, for finite clusters the boundary conditions are crucial. For the chains 
studied here when N = 4n + 2 (where n is a natural number) the periodic boundary 
conditions (PBC) are used, whereas for N = An the antiperiodic boundary conditions 
(ABC) lead to a lower-energy state. Namely, the terminal-atomic annihilation operators 
should be defined as follows 

Ojv+v — > a la forPBC, 
ajv+i,cr _a io for ABC. 

Therefore, the terminal hopping term involving the end atoms changes sign for ABC, 
while the interaction terms ~ n^riii, UiUj and Sj ■ Sj remain unaltered. Obviously, the 
periodic boundary conditions are proper for any N, when we have a ring geometry. The 
situation becomes more involved when N is an odd number. Namely, we write down 
the Hamiltonian (jfiHj) in the form 

H = e e J f J2 n ^ + t J2 e~ im a] a a j+1(J + Uj^ n^n u + £ K^SriMj, (78) 

i(T jtr i i<j 

where $ is the fictious (dimensionless) flux through the ring. One can show that the 
unitary transformation Cj a — ► e~^l^ Cj„ allows for accumulation of complex phase factor 
in the terminal hopping term, which then takes the form t{eT % ® a} N a a\ a + h.c.) and 
this can be regarded as generalized boundary condition a^+iu — » e l ®' N ax a . We do 
not distinguish between the system with a fictious flux and the generalized boundary 
conditions. The presence of the flux can be regarded as an accumulation of the Berry 
phase during the motion of individual electron in the milieu of all other electrons. With 
such interpretation the boundary conditions apply also to linear chain. 

The electron momentum for nanochains with such boundary conditions is displayed 
in Fig. ESI The discrete momenta, corresponding to the solution of the single particle 
part of (|7Hj) for a finite N, are given by 

2nn - $ . , 

K= N ,0<n<N (79) 

The optimal BCs, corresponding to the minimal ground-state energy Eg, with respect 
to $, are realized for $ = when iV = 6,10,14,... (periodic BCs), $ = 7r when 
N — 4, 8, 12, . . . (antiperiodic BCs), and $ = 7r/2, 3/2tt when iV is odd. A basic analysis 
of Eq. (fTSj) shows, that for the optimal BCs, the Fermi momentum value kp = ir/(2R) 
is never reached for even N, whereas for odd N it happens for a single value of n. This 
circumstance has important implication for the nanochain electronic structure, however, 
almost does not affect its transport properties, as discussed below. The exact meaning 
of the averaging procedure when drawing Fig. EHl is elaborated elsewhere. Analogously, 
the effect of BCs on the spin-spin correlation function < Sj • Sj > is also touched upon 
there. 
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An explicit form of the dispersion relation derived from the spectral density 
functions with inclusion of the generalized boundary conditions, is shown in Fig. |22J 
We see again, that the spin splitting is present. Also, one can see the difference between 
either dispersion relation for noninteracting particles or with that calculated in the 
Hartree-Fock approximation for the Slater antiferromagnetic chain. 



5. Transport in nanosystems 

In this Section we complement the discussion of nanochains by providing the values of 
Drude weight (charge stiffness) and optical gap. We also show, how the EDABI method 
can be used to calculate the conductance of an open system: a quantum dot attached 
to the leads. 



5.1. Drude weight for nanochains 

The real part of the optical conductivity at zero temperature is determined by the 
corresponding real part of the linear response to the applied electric field |3*T| . and can 
be written as cr(u) = D5{uj) + cr rog (c<;), where the regular part is 



= I E K *; ljpl ! 0)l S (<" - ( E n -Eo)), (80) 
N 5£ E n - Eo 

whereas the Drude weight (the charge stiffness) D is defined by 

D = -— <* | T |* > ~-Y (81) 
AM 01 1 0/ N E n -E ' 1 ; 

where the kinetic-energy term T is the same as the second term in Eq. 1)780. and 
the diamagnetic current operator is defined as usual: j p = itJ2j a ( a ja a j+^ — h.c). 
Here the states \^ n ) in Eqs. (jHU|) and (jHIJ) are the eigenstates of the Hamiltonian (f7Hj) 
corresponding to the eigenenergies E n , again with boundary conditions which minimize 
the ground-state energy for a given system size N (see previous Section). We calculate 
matrix elements (\& n | j p \ ^f ) numerically by adapting the method developed by Dagotto 
[32] . which is stable for even number of electrons only. Namely, we diagonalize the 
Hamiltonian (|78j) in the Fock space with the help of Davidson technique |HH! , and then 
repeat the procedure starting from a specially prepared initial state j p \ ^o). Replacing 
the Lanczos scheme, utilized in Ref. [32J, by the Davidson one provides an excellent 
accuracy for either even or odd number of electrons. The evolution of Drude weight 
with iV and R is shown in Fig. EH In the half-filled case N e = N (cf. Fig. 123k ) 
Drude weight gradually decrease with N, as we have shown for even N |34j . The most 
interesting feature of these results is, that the curves for odd N fits smoothly between 
those for even N, with very weak parity effect, totally incomparable with that present 
in the charge gap and electron momentum distribution [33]. This observation can be 
understood when we take into account, that the Drude weight defined by Eq. ED is the 
integral quantity, involving the summation over all the excited states of the Hamiltonian 
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Figure 23. Electron momentum distribution for chains of N = 8 -j- 13 atoms: (a) 
datapoints for optimal boundary conditions (BCs) and the sample n(k) curves averaged 
over BCs (solid lines); (b) the original n(k q (4>)) functions used for the averaging. 



Eq. ()78|) . so it cannot be determined only by the electronic structure near the Fermi 
points, particularly for a small system. 

The parity effect on Drude weight disappears for the system with a single hole 
(N e = N — 1, cf. Fig. l2"4"b). in which the magnetic frustration is absent. For this 
case, the Drude weight evolution with R is very interesting. In the weak-correlation 
range (R/ao < 2) the chain shows a highly-conducting behavior for each N. Next, in the 
intermediate range (R/a Q = 4 -=-5) the Drude weight decrease rapidly with N, indicating 
an insulating (Mott-Hubbard) state in the large- N limit. In the strongly-correlated 
range (R/ao ~ 10) the Drude weight approaches again its maximal value D — 1. Such 
a behavior can be explained, when we analyze the situation in two steps: First, for low 
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INTERATOMIC DISTANCE, R (a.u.) 

Figure 24. Normalized Drude weight for nanochains in the half-filled case (a), and 
for a system with a single hole (b) . The optimal boundary conditions are specified for 
each dataset. 

values of R, the bandwidth-to-interaction ratio is small, and the system with a single 
hole does not differ significantly from a half-filled one. This is why in both cases Drude 
weight decreases gradually with both N and R, as the tunneling amplitude through the 
barrier of a finite width. Second, for the largest values of R, the system can be described 
by an effective t — J model [21] with a coupling constant J = 4t 2 /(U — K) <C \t\ (where 
K = Kjj + i denotes the nearest neighbor Coulomb repulsion), which corresponds to an 
asymptotically-free hole motion. Then, it become clear that in the intermediate range 
the Drude weight has to suppressed, what can be interpreted in terms of a partially 
localized spin-ordered state. It would be very interesting to test experimentally this 
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Figure 25. Optical gap for nanochains in the half-filled case. 



result, possibly for a mesoscopic atomic ring. 

The smooth behavior of the system transport properties with N is illustrated in Fig. 
1251 on the example of other relevant quantity, the optical gap AE opt = E^ — Eq (where E^ 
is the energy lowest-lying excited level for which (^i\j p \^/ ) ^ 0). Again, adding a 
single hole totally removes the weak parity effect (not shown) . The large- R behavior of 
AE opt could be explained by the value for an insulating limit AE opt « U — K, whereas 
small R we have a strong N dependence similar to those observed for D. The later 
shows, that the contribution of the state to the sum (}8l"|) is dominating for such a 
partly-localized quantum nanoliquid 

5.2. H 2 molecule as a quantum dot 

We have recently shown, that EDABI method is useful to investigate the zero- 
temperature conductance of diatomic molecule, modeled as a correlated double quantum 
dot attached to noninteracting leads jl]. This became possible when utilizing the Rejec- 
Ramsak formulas, relating the linear-response conductance to the ground-state energy 
dependence on magnetic flux j2Zj. Here we complement the discussion with the stability 
analysis showing, that large coupling between molecule and the leads may provide the 
possibility for interatomic distance manipulation in an experiment. 

We start with the Hamiltonian which is a generalization of the Anderson impurity 
model [HB], and can be written as 

H = H L + V L + H C + V R + H R , (82) 

where He models the central region, Hl(r) describes the left (right) lead, and Vl{r) is 
the coupling between the lead and the central region. Both Hur) and Vl(r) terms have 
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Figure 26. Diatomic molecule as a double quantum dot. 



a tight-binding form, with the hopping t and the tunneling amplitude V 

H L (R)=t ^2 ( a ]* a j+i,v + n - c -)> (83) 

V L (R) = V ^2(a[ {R)a a 1{2)a + h.c.), (84) 

17 

as depicted schematically in Fig. [2I3 The index j = L(R) corresponds to the last atom 
of the left (right) lead. The central-region term He is a version of Hamiltonian ("o"8~|) for 
two atoms, and describes a double quantum dot with electron-electron interaction 

H c = (e a - eV G ) ^ n i ~ f ' ^ ( a L a 2a + n -c) 

j=l,2 <j=t,J. 

+ C/ ^2 n i\ n H + ^1^2 + (Ze) 2 /R, (85) 

where e a is atomic energy, Vq is an external gate voltage, t' is the internal hopping 
integral, U and K represents the intra- and inter-site Coulomb interactions, respectively, 
and the last term describes the Coulomb repulsion of the two ions at the distance R. 
Here we put Z = 1 and calculate all the parameters e a , t', U, and K as the Slater 
integrals for ls-like hydrogenic orbitals \?i a (r) = exp(— a|r|), where a -1 is 

the orbital size (cf. Fig. I2fij) . The parameter a is optimized to get a minimal ground 
state energy for whole the system described by the Hamiltonian ([82)1 . Thus, following 
the idea of EDABI method |H4j . we reduce the number of physical parameters of the 
problem to just a three: the interatomic distance R, the gate voltage Vq, and the lead- 
molecule coupling V (we put the lead hopping t — 1 Ry = 13.6 eV to work in the 
wide-bandwidth limit). 

We discuss now the molecule conductivity calculated from the Rejec-Ramsak two- 
point formula [3*7] : 

G = G'osin 2 7r[ J E;(7r)- J E(0)]/2A, (86) 

where Go = 2e 2 /h is the conductance quantum, A = 1/Np(ep) is the average level 
spacing at Fermi energy, determined by the density of states in an infinite lead 
p(ep), E(tt) and E(0) are the ground-state energies of the system with periodic and 
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Table 11. Same as in Tabic calculated in Gaussian-type STO-3G basis, with the 
inclusion of the long-range Coulomb interactions. 
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Figure 27. Zero-temperature conductance of the system in Fig. |2EI as a function 
of the gate voltage Vq and interatomic distance R. Is orbital size a -1 is optimized 
variationally. The lead parameters are t = lRy and V — O.bt. 
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Figure 28. The optimized inverse inverse orbital size (a) and average central region 
occupancy (b) for the diatomic molecule attached to noninteracting leads characterized 
by the parameters t — lRy and V = 0.5i. 



antiperiodic boundary conditions, respectively. E{4>) is calculated for = 0, n within 
the Rejec-Ramsak variational method [HI], complemented by the inverse orbital size a 
optimization, as mentioned above. We use typically iV = 10 2 -j- 10 3 sites to reach the 
convergence. 

In Fig. E3 we show the conductivity for V = 0.5t, and different values of the 
interatomic distance R. The conductance spectrum evolves from the situation well- 
separated peaks corresponding to the independent filling of bonding and antibonding 
molecular orbitals (R < 2a , where a is the Bohr radius), to the single peak in the 
intermediate range (R ~ 3ao), which decays when t' < K for large R. 

Probably, the most interesting feature of the conductance depicted in Fig. Ellis their 
strong asymmetry for small R. Namely, the low-V^ conductance peak, corresponding 
to the system filling (711 + 712) ~ 1 (one hole) is significantly wider than the high- Vq 
peak for (ni+n 2 ) ~ 3 (one extra electron). Such a particle-hole symmetry breaking is 
the novel feature, observed when including the correlation-induced basis optimization, 
and absent in the parametrized-model approach [3T[ I30*j . It is also a new feature of 
a nanosystem, do not observed in mesoscopic double quantum dots [41J, where the 
particle-hole symmetry holds. 

The relation between the observed asymmetry, basis renormalization, and electron 
correlations can be clarified as follows. First, one can observe that the optimal values 
of the variational parameter a, provided in Fig. l2"Hk. decrease dramatically for high Vq, 
corresponding to the overdoped situation (rii + n 2 ) > 2 (c/. Fig. I2%b). This is because 
the system minimize energy of double occupancies, which is of the order U ~ a [39|. 
Then, we focus on an small R limit, in which the well separation of molecular orbitals 
allows one to approximate the expression for a single impurity at zero temperature [12] 

G = G sin 2 (n(n k }/2), (87) 

where (n^) is the bonding (k = 0) or antibonding (k = it) orbital occupancy. Expanding 
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Table 12. The binding energy AE and the bond length i? m i n for different coupling 
to the lead V . The conductance at the energy minimum is also provided. 



v/t 


AE 


-Rmin/ 


^min /^O 


0.0 


-0.296 


1.43 





0.1 


-0.293 


1.44 


10~ 4 


0.2 


-0.282 


1.46 


0.004 


0.3 


-0.180 


1.52 


0.023 


0.4 


-0.040 


1.59 


0.093 


0.5 


+0.074 


1.93 


0.422 



around the maximum we have G(Vg) — C(V^' fc ) ~ —G (x c ' k /^) 2 (Vg — V^*), where V G ,k 
(k = 0, 7r) is the low/high voltage peak position and Xc is the charge susceptibility, 
which may be approximated as Xc ~ d(n\ + n2) /9Vg, since the orbitals are filled 
separately. Values of the later derivative read from Fig. l2*Hb around the low and high 
voltage peak positions ((ni+na) ~ 1 and w 3, respectively) provides a clear asymmetry. 
Moreover, the expansion of G(Vg) allows one to roughly estimate the peak width as 
AVg ~ Xc 1 ~ {U + K)/2 ~ a, what provides an indirect correspondence between the 
spectrum asymmetry and basis renormalization. 

The practical possibility of the conductance measurement involving atoms 
manipulation (changing of R) were also explored in terms of the system stability. 
Namely, the binding energies AE = E Rinin — E^ (where i? m i n is the bond length), 
listed in Table El shows that the system become metastable around V = 0.5t (AE > 
indicates a local energy minimum). Therefore, the binding of the atoms to the lead 
become stronger than interaction between the atoms, that may allow individual atoms 
manipulation. 

We would like to stress here, that ab initio analysis have never been performed for 
the open system with strong electron correlation. Apart from simplifying the discussion 
(the model parameters are calculated as a function of interatomic distance), such an 
approach leads to new physical effects, since the particle-hole symmetry is broken, as 
exemplified on the example of the conductivity calculations. 

6. Concluding remarks 

In this paper we have provided a detailed discussion of the method based on the 
diagonalization in the Fock space, when combined with with an ab initio adjustment of 
single-particle wave functions in the interacting state of N-particle system (the EDABI 
method). The method can be improved upon by systematically increasing the number of 
wave functions in the basis [wi(r)] when defining the field operator for the system (thus 
enriching the model). We have illustrated our method with the discussion of several 
nanoscopic systems ranging from atoms to nanochains and quantum dots. 
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The method is useful when the exact solution of a model is available in the Fock 
space. Such situation happens also for one dimensional atomic chain represented by 
the Lieb-Wu solution of the Hubbard model [30! • I n that situation, exact Wannier 
or Bloch wave functions can be constructed. Other applications such as solving the 
magnetic impurity in a nonmetallic environment are also possible, although not carried 
out explicitly as yet. 

The importance of our approach is, among others, in showing explicitly that the 
concept of the statistical distribution of particles with quasimomentum, as a good 
quantum number, is feasible for a relatively small interatomic distances. This means, 
that the corresponding N-electron states form a quantum nanoliquid for not too-large 
inter-distance in nanowire. In connection with this, in replying the question how small 
a piece of metal can be, we can say that the nanoliquid exhibits metallic conductivity 
for eV > AE(N,R), where V is the voltage applied to the system and AE(N,R) is 
the energy difference between the highest occupied and the lowest unoccupied state, i.e. 
that the momentum is a good quantum number then. For an intermediate interatomic 
distace, we have a gradual transition to a Mott insulator, above which the monoatomic 
(quantum) nanowire is useless for electronic applications. 
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